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MODELING OF UNIT-CELLS WITH Z-PINS USING FLASH: 
PRE-PROCESSING AND POST-PROCESSING 

Ronald Krueger 1 


ABSTRACT 

Although the toughening properties of stitches, z-pins and similar structures have 
been studied extensively, investigations on the effect of z-pins on the in-plane 
properties of laminates are limited. A brief summary on the effect of z-pins on the 
in-plane tensile and compressive properties of composite laminates is presented 
together with a concise introduction into the finite element code FLASH. The 
remainder of the report illustrates the modeling aspect of unit cells with z-pins in 
FLASH and focuses on input and output data as well as post-processing of results. 
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Yp 

shear strain 

Yy 

Yield strain in shear 

A 

loading parameter 
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Lamina Poisson’s ratio 

4> 

Fiber misaligmnent angle 
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Strength, critical value of stress 
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Yield strength in tension 
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Compression strength of skin/ stiffener- flange laminate 
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Normal stress in 1, 2 direction 

a 12> a 21 

Shear stress in 1-2 plane 

X e 

Effective yield strength in shear 

Xy 

Yield strength in shear 

T X y, T'yx 

Shear stress in X-Y plane 

t 12’ r 21 

Shear stress in 1-2 plane 


1. BACKGROUND 

One of the most common failure modes for composite structures is delamination [1-3]. 
The remote loadings applied to composite components are typically resolved into interlaminar 
tension and shear stresses at discontinuities that create delaminations under mixed-mode 
conditions. In the past research has focussed on stitching to increase the delamination resistance 
of composites [4-8]. More recently z-pins 2 have been proposed to stitch the material together 
through a combination of friction and adhesion [9-12]. Z-pins are pultruded rods of carbon fiber 
and epoxy matrix with an approximate diameter of 280 pm or 508 pm. The z-pins are inserted 
through the thickness of a laminated composite, which is then autoclaved to cure the epoxy as 
normal. This approach to through-thickness reinforcement offers an alternative to stitching, and 
can provide much higher areal densities of reinforcement [13]. 

Examination of the literature shows that the toughening properties of stitches, z-pins and 
similar structures have been studied extensively, however, investigations on the effect of z-pins on 
the in-plane properties of laminates are limited [13-17], Steeves examined the effect of z-pins on the 
in-plane tensile and compressive properties of composite laminates [13]. Disruption in the 
alignment of the fibers in a fiber composite leads to a significant reduction in the in-plane 
compressive strength [18]. Since the diameter of the z-pins (~280 pm) is large relative to the 
diameter of the fibers (~7 pm) in the composite, the z-pins may cause significant misalignment of 
the fibers of the composite. A sketch of the typical distribution of fiber misalignment around a z-pin 
is shown in Figure 1. The z-pin is surrounded by a resin-rich pocket on two sides, caused by the 
fibers deflecting around the z-pin and leaving a gap in the composite. The maximum misalignment 
of the deflecting fibers occurs in four locations on the four flanks of the resin-rich pocket. 


2 The generic term z-pin will be used throughout the paper versus the trade mark Z-Fiber™ registered by Aztex Inc. 
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Steeves’ experiments and numerical finite element simulations using FLASH showed that 
z-pins reduce the in-plane strength of composite laminates. His experimental data from 
compressive tests showed that the insertion of a field of z-pins into a laminated fiber composite 
reduces the in-plane compressive strength of the composite by 30% or more. He assumed that 
the reduction in compressive strength is a consequence of the increased fiber waviness caused by 
the z-pins [13]. Compression tests within the scanning electron microscope showed that the z- 
pinned specimens develop microbuckles in the region of greatest fiber waviness around the resin- 
rich pocket, as would be predicted by the microbuckling model of Budiansky and Fleck [19]. 

The objective of this report is to present a brief introduction into the finite element code 
FLASH and illustrate the modeling aspect of unit cells with z-pins in FLASH. The report focuses 
on the detailed description of input data for FLASH and output obtained from the analysis as 
well as the post-processing of results. The report is intended to supplement existing user 
manuals [20, 21], First, input is described for a finite element model and material data from 
reference [13]. Second, the correct load input and boundary conditions for compression, tension, 
shear and combined loading is illustrated. Third, for a composite made of carbon/epoxy modeling 
for large (0.508 mm diameter) and small (0.28 mm diamter) z-pins is described for 2% and 4% z- 
pin areal density. Input data are determined for unit cells subjected to five different load cases 
such as axial compression, combined axial compression and transverse tension, combined axial 
compression and shear loading as well as combined axial compression, transverse tension and 
shear loading. Additionally, the output obtained from FLASH is explained and separate post- 
processing options are discussed. The interpretation of the output with respect to the strength of 
z-pin reinforced lamina is discussed in detail in a separate report [22], 


2. INTRODUCTION TO FLASH 

In order to better assess the influence of critical parameters on lamina compression 
strength. Fleck and Shu developed a finite element code called FLASH [20, 21, 23]. The FE code 
is based on a 2D general Cosserat couple stress theory that assumes the unidirectional composite 
lamina is a homogeneous anisotropic material that carries couple stress as well as classical 
Cauchy point stress [24, 25]. The constitutive response is deduced from a unit cell consisting of 
a fiber, represented by a linear elastic Timoshenko beam, embedded in a non-linear elastic-plastic 
matrix [26]. The fiber diameter, d, is the length scale in the constitutive law that controls fiber 
bending resistance. The continuum theory was implemented within a two-dimensional finite 
element code that uses 6-noded triangular elements with 3 degrees of freedom at each node (two- 
displacements and one rotation corresponding to rotation of the fiber cross section). The finite 
element procedure is based upon a Lagrangian formulation of the finite deformation of the 
composite and can accommodate both geometric and material non-linearities [23]. The code 
models finite deformation using a Newton-Raphson incremental solution procedure with a 
modified Riks algorithm in the final stage to handle snap-back behavior associated with fiber 
micro-buckling [27]. Boundary loading is piecewise proportional with a loading parameter, A, for 
each loading stage [20, 21, 23]. 

The FLASH code assumes micro-buckling initiates from an imperfection in the form of 
fiber waviness. Inputs include material properties (volume fraction V f , and stiffness properties 
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E l , Ej, Ej c , G[ [, Gf, all normalized by the shear yield strength, and the Ramberg-Osgood 
strain hardening law parameters (a, n). FLASH allows options for input of fiber misalignment 
angle due to fiber waviness either as (1) an elliptical patch of waviness, or (2) an arbitrary 
distribution of initial fiber waviness through initial misalignment angle, (p, at the Gauss 
integration point for each element [20, 21, 23]. The first option prescribes the elliptical patch 
along one edge of the unit cell. Steeves used the second option to input fiber misalignment 
distribution obtained from analysis of digital images of specimens taken in a Scanning Electron 
Microscope (SEM) [13]. 


3. GENERIC PRE-PROCESSING 
3.1 Unit Cell Geometry 

The z-pin reinforced composite is divided into unit cells as shown in Figure 2. The size of 
the unit cell depends on the areal density, r z (in %) of the z-pins and the diameter of a single z- 
pin, D z as shown in Figure 2. The spacing L Z =H Z for a perfect, rectangular z-pin field can be 
calculated as 



where A z is the cross sectional area of a single z-pin 



and A r denotes the fraction of the total reinforced area covered by z-pins 



The length of the resin pocket, C, may be determined from micrographs of the reinforced laminate 
as shown in Figure 3. 


3.2 Input of Geometrical Points 

For the analysis with FLASH only one unit cell is modeled as shown in Figure 4. The resin-rich 
pocket and the z-pin are modeled as voids within the material. Considering axial compression (x- 
direction) as the driving failure load, modeling a void is justified because the matrix strength and 
stiffness are low relative to the strength and stiffness of the composite, and the resin-rich pocket 
can considered to be significantly cracked early in the simulation. Hence, the simulation models 
the unit cell of a composite including fiber misalignment and internal holes [13]. All meshes 
generated were composed of six-noded triangular plane-strain elements. The size of the elements 
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was varied to provide the greatest mesh refinement near the resin pocket, and in the region of 
greatest fiber misalignment. 

First geometric parameters were taken from reference [13]. Since the fiber is completely 
surrounded by resin the transverse dimension of the modeled void, D ’ z , is increased by 0.02 mm 
compared to the z-pin diameter, D z as shown in Figure 3a. The finite element code FLASH 
requires the input of geometric data normalized with the fiber diameter of the composite, d as 
listed in Table 1. Geometrical data for the areas to be meshed is input using the *geometrical 
points command in FLASH. The coordinates of the points used to create the finite element 
model of the unit cell with z-pin are printed below and their physical location is sketched in 
Figure 5. 


* GEOMETRICAL POINTS 
26 


1 

0.00 

0.00 

2 

125.00 

0.00 

3 

250.00 

0.00 

4 

0.00 

103.60 

5 

125.00 

103.60 

6 

250.00 

103.60 

7 

110.00 

105.00 

8 

140.00 

105.00 

9 

0.00 

124 . 00 

10 

33.00 

124 . 00 

11 

217 . 00 

124 . 00 

12 

250.00 

124 . 00 

13 

32.00 

125.00 

14 

217 . 99 

125.00 

15 

0.00 

126.00 

16 

33.00 

126.00 

17 

217.00 

126.00 

18 

250.00 

126.00 

19 

110.00 

145.00 

20 

140.00 

145.00 

21 

0.00 

146.40 

22 

125.00 

146.40 

23 

250.00 

146.40 

24 

0.00 

250.00 

25 

125.00 

250.00 

26 

250.00 

250.00 


The geometric parameters were used to generate the finite element mesh of the unit cell including 
the resin-rich pocket as shown in Figure 4. 
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3.3 Input of Mesh Data 

Mesh seed data for the areas to be meshed is input using the *mesh sets command in 
FLASH. The input to create 10 mesh sets which form the finite element model of the unit cell 
with z-pin (see Figure 4) are printed below. The physical location of the mesh sets is sketched in 
Figure 6. 


*MESH SETS 

10,6 

1 , 1 , 4 , 2 , 2 , 5 , 8 , 10 , 0 . 1 , 0 . 1 , 0 . 1 , 0.1 

2 , 2 , 5 , 3 , 3 , 6 , 8 , 10 , 0 . 1 , 0 . 1 , 10 , 10 

4 , 4 , 9 , 5 , 7 , 10 , 8 , 10 , . 25 , . 25 , 0 . 1 , 0.1 

5 , 8 , 11 , 6 , 6 , 12 , 8 , 10 , . 25 , . 25 , 10 , 10 

9 , 9 , 15 , 10 , 13 , 16 , 3 , 10 , 1 , 1 , 0 . 1 , 0.1 


14 , 

17 , 

12 , 

12 , 

18 , 

3 , 

10 , 

1 , 

1 , 

10 , 

10 


15 , 

21 , 

16 , 

19 , 

22 , 

8 , 

10, 

4 , 

4 , 

0 . 1 , 

0 

. l 

20 , 

22 , 

18 , 

18 , 

23 , 

8 , 

10, 

4 , 

4 , 

10 , 

10 


21 , 

24 , 

22 , 

22 , 

25 , 

8 , 

10, 

10 , 

, 10 

, 0. 

1 , 

0.1 

22 , 

25 , 

23 , 

23 , 

26 , 

8 , 

10, 

10 , 

, 10 

, 10 

r 

10 


3.4 Input of Misalignment Data 

The finite element code FLASH also calls for the input of the fiber misalignment angle [20, 
21], The misalignment angle is input using the ^misalignment type and the *misalignment 
command in FLASH. The input used for all analysis is printed below. 

*MISALIGNMENT TYPE 
2 

*MISALIGNMENT 

74 , 74 , 3 . 425 , 3.425 

Preliminary analysis with N=\ as input for *misalignment type which was supposed 
to generate an elliptical patch of waviness was not successful. Post-processing - discussed later - 
only showed a small band of waviness along the left edge of the specimen (x=0). The waviness in 
the rest of the specimen was zero. Therefore this approach was abandoned and JV=2 was chosen 
which requires that the discretized misalignment data is provided in a separate input file 
(phiex.dat). Input of an arbitrary distribution of the fiber misalignment is possible as by 
sampling a micrograph as shown in Figure 7 [13]. However, if this information is not available, 
misalignment angles may be chosen to study the misalignment effect on the reduction of in-plane 
strength. 

To illustrate the input an example is shown in Figure 8 for 10 by 10 sample points on the 
unit cell which corresponds to a matrix of 10 by 10 data points in the phiex.dat file. The 
distance of one point to the adjacent points was set to Ax=26.0 units in x and Ay= 26.0 units in y- 
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direction which extends beyond the unit cell (250 by 250 units) modeled. The distance of points 
has to be chosen so that the field of data points is larger than the unit cell. The user is advised to 
check the misalignment input after the execution of prep has generated the flash . inp file and 
adjust the distance of adjacent points accordingly. If the misalignment data originates from 
sampling a micrograph the sample area has to be larger than the model. See reference [21] for 
details. Only the file named phiex . dat will be used during the initial step of the analysis when 
the file flash . inp. is generated during the execution of prep . 

For the current analyses constant fiber misalignment values of 0°, 1°, 2°, 3°, 4°, 5°, 6°, 
7°, 8°, 9°, and 10°, were chosen and 11 phiex . dat files were created. Each phiex.dat file 
contains a field of 74 by 74 data points which corresponds to a matrix of 74 by 74 sample points 
on the unit cell as shown for the example of 5° fiber misalignment. 
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5.000 

5.000 

5.000 

5.000 

5.000 

5.000 








The distance of one point to the adjacent points was set to 3.425 units in x and 3.425 units in y- 
direction to cover the entire unit cell (250 by 250 units). 


3.5. Input of Constraints 

As mentioned above the composite is divided into unit cells, and only one cell is modeled 
using periodic boundary conditions. For these models, periodic boundary conditions are 
prescribed by applying compressive loading symmetrically to the specimen, and constraining 
three degrees of freedom to prevent rigid body translation and rotation as shown in Figure 9. The 
constraints are input using the *type 1 constraints command in FLASH. The input used for 
all analysis is printed below. 

*TYPE 1 CONSTRAINTS 
3 

1 , 1 , 2 , 21 , 1 , 1 , 0 
9 , 9 , 15 , 4 , 4 , 2 , 0 
24 , 24 , 25 , 21 , 1 , 1 , 0 
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3.6. Input of Load Data 

The loading was input using the *stress loading command in FLASH. The axial 
compression stress sketched in Figure 9 is gradually incremented by FLASH until it reaches the 
specified limit defined by the user ( o xx /r y =-1000). The limit was deliberately chosen to be well 
above failure, i.e. formation of kink bands and fiber microbuckling to assure that the analysis 
reached the failure point and did not terminate at a chosen lower stress [20]. The input used for 
axial compression loading is printed below. 


*STRESS LOADING 


10 , 1 


1 , 1 , 

4 , 1 , 1 

- 1000 , 

0 , 0 , 0 , 

4 , 4 , 

9 , 1 , 1 

- 1000 , 

0 , 0 , 0 , 

9 , 9 , 

15 , 1 , 1 

- 1000 , 

0 , 0 , 0 , 

15 , 15 

, 21 , 1 , 

- 1000 , 

0 , 0 , 0 , 

21 , 21 

, 24 , 1 , 

- 1000 , 

0 , 0 , 0 , 

3 , 3 , 

6 , 1 , 1 

- 1000 , 

0 , 0 , 0 , 

6 , 6 , 

12 , 1 , 1 

- 1000 , 

0 , 0 , 0 , 

12 , 12 

, 18 , 1 , 

- 1000 , 

0 , 0 , 0 , 

18 , 18 

, 23 , 1 , 

- 1000 , 

0 , 0 , 0 , 

23 , 23 

, 26 , 1 , 

- 1000 , 

0 , 0 , 0 , 


0 

0 

0 

1 

0 

1 

0 

0 

0 

1 

0 

1 

0 

1 

0 


3.7. Requests for Output 

During the analysis stresses, strains and total fiber alignment angle is written to the file 
flash. out after every N load steps. The data output request is managed using the *output 
frequency commands in FLASH as printed below for an output request every 1000 steps. 

* OUT PUT FREQUENCY 
1000 

During the analysis, for three boundaries (as shown in Figure 9) data are output to the 
output file post . d after each load step. For boundary sections 1 and 2 the net force on each 
boundary section, resolved along the nominal fiber direction, are given as output. Boundary 3 is a 



single node and the displacements of this node are given as output. The data output request is 
managed using the *output list command in FLASH as printed below. The boundary sections 
selected are shown in Figure 5. 

* OUT PUT LIST 
3 

1 , 1 , 4 , 1 , 1 

1 , 1 , 4 , 1 , 1 

1 , 1 , 4 , 1 , 1 

Additionally, the value of the global loading parameter A, is written to the output fde at 
each loading step. 

The format of the output files and post-processing of the results is discussed in a 
separate section below. 


3.8. Input of Material Data 

The finite element code FLASH requires the input of material data normalized with the 
composite yield shear strength r y as given in Table 2 for a generic graphite epoxy material [13]. 
Additional material parameters used by FLASH are the constants a=0.429 and n = 3 in the 
Ramberg-Osgood strain hardening law, and the fiber volume fraction F/ =0.63 of the composite. 
The material data is input using the *material constants and *vfraction commands in 
FLASH as printed below. 

* MATERIAL CONSTANTS 
7 

1200 , 175 , 74 , 41 , 0 . 429 , 3 . 0 , 87 
*VFRACTION 
0.63 


3.9. Input for a Finite Element Model without Z-Pin 

For reference purposes a unit cell without z-pin was modeled and analyzed with FLASH. 
The mesh with applied load and boundary condition is shown in Figure 10. The geometrical points 
and the required input data to generate the mesh have been included in Figure 10. 


4. PREPROCESSING FOR COMBINED LOADING 

For the initial analyses performed, load and boundary conditions were taken from 
reference [13]. Compressive loading was applied symmetrically to the specimen, and three 
degrees of freedom were constrained to prevent rigid body translation and rotation as shown in 
Figure 9. 

In general, z-pin reinforced composite structures are subjected to multi axial loads and 
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input stresses in the unidirectional ply modeled in FLASH need to be calculated. For composite 
plates subjected to external compression and combined compression and shear loading input 
stresses for FLASH were calculated. Using classical laminate theory the external loads were 
resolved into stresses for the individual plies as shown in Figure 1 1 . Calculated stresses in the 
critical unidirectional ply consist of axial compression, transverse compression and tension as well 
as shear stresses. For these combined loading conditions the above mentioned boundary 
conditions and load input had to be modified. 


4.1. FE Models of Z-Pin Unit Cells Subjected to Axial and Transverse Loading 

The finite element mesh with load and boundary conditions for axial and transverse 
compression is shown in Figure 12 with the corresponding input data. Based on results from 
classical laminate theory mentioned above, the transverse compression load was chosen to be 2% of 
the compression load applied axially. In order to avoid the collapse of the void under transverse 
compression points A and B on the periphery of the void could be constrained in transverse 
direction as shown in Figure 12. In reality, the boundary condition of zero displacements on the 
periphery of the void is not met due to pin compliance. However, the assumed boundary 
condition is adopted in order to obtain an estimate for the extreme case of a rigid pin [28], 
Nevertheless, analysis without constraints suggested that constraining points A and B for small 
transverse compression loads leads to overly conservative results. 

The finite element mesh with load and boundary conditions for axial compression and 
transverse tension are shown in Figure 13. As for the previous case the transverse load was chosen 
to be 2% of the compression load applied axially. For transverse tension it is not required to avoid 
the collapse of the void. 


4.2. FE Models of Z-Pin Unit Cells Subjected to Axial Compression and Shear Loading 

The finite element mesh with load and boundary conditions for axial compression 
combined with shear are shown in Figure 14. Based on results from classical laminate theory 
mentioned above, the shear load was chosen to be 10% of the compression load applied axially. For 
the combined axial compression and shear load case it was decided not to prevent the collapse of 
the void as explained in detail in the appendix of reference 23. 

The finite element mesh with load and boundary conditions for axial compression and 
transverse tension combined with shear are shown in Figure 15. As for the previous case the shear 
load was chosen to be 10% of the compression load applied axially and the transverse tension load 
corresponded to 2% of the compression load as discussed in the previous section. For this 
combined load case it was decided not to prevent the collapse of the void as explained in detail in 
the appendix of reference 23. 


5. UNIT CELL MODELS FOR 2% AND 5% Z-PIN REINFORCEMENT. 

For another carbon fiber reinforced epoxy material the data used as input for the analyses 
with FLASH are listed in Table 3. The specimens were reinforced with small z-pins of 
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D z =0.279 mm (0.011 in) and large z-pins of D z =0.508 mm (0.02 in) diameter. Three specimen 
types were manufactured containing reinforcement fields with r z = 2% areal density for the large 
z-pin and r 7 = 2% and r z =4% areal density for the small z-pins respectively. Geometric data such 
as z-pin spacing, L Z =H Z , as well as width, D ’ z , and length, C, of the resin pocket as defined in 
Figure 3 was measured from micrographs taken of different specimens. The averaged data used as 
input for the analysis is summarized in Table 4. 

Based on results from classical laminate theory shown in Figure 11 and Table 5, five load 
cases were considered as input for analysis of the unit cells: A pure axial compression load case, a 
combined axial compression - 2% transverse tension load case, a combined axial compression - 
10% shear load case, a combined axial compression - 50% shear load case, and an axial 
compression - 2% transverse tension load case combined with 10% shear loading. The input for all 
five cases is summarized in Table 6. 

The finite element meshes are shown in Figures 16, 17, and 18 for the small pin with 2% 
and 4% areal density and the large pin with 2% areal density. As before the resin-rich pocket and 
the z-pin are modeled as voids within the material. It is assumed that the fiber is completely 
surrounded by resin as shown in Figure 3 and therefore the transverse dimension of the modeled 
void, D 2, is increased by 0.02 mm compared to the z-pin diameter, D z . The normalized dimensions 
required as input were included in the figures as well as the input for the *geometrical points 
command in FLASH. Loads given in Table 6 were applied as discussed in section 4. 


6. POSTPROCESSING OF RESULTS 
6.1. FLASH Output Files 

During the analysis nodal point coordinates, element topology and results are is written 
to the file flash . out after every N load steps. The file flash . out contains a header, 

lCcomnew 0.00000E+00 

the nodal point coordinates 
2C 

-1 1 0 . 00000E+00 

-1 2 0 . 00000E+00 

-1 3 0 . 00000E+00 

-1 4 0 . 00000E+00 

-1 ... 

-3 

the element topology 
3C 

-118 

-2 1 374 372 155 373 154 

11 


FIBER MICROBUCKLING 3 


0 . 00000E+00 
0 . 15646E+02 
0 . 31293E+02 
0 . 42553E+02 



374 


2 156 155 


-12 8 
-2 13 

-13 8 

-2 372 710 708 581 709 580 

-1 ... 

-3 

the initial fiber misalignment (f> 

-3 

100CLCASE1 O.OOOOOE+OO Fiber Misalignment 3 0 


-4 

MISALIGN 1 

3 

0 

-5 

PHIBAR 1 

1 

1 

-1 

1 

8 0 

0 

6 

-2 

1 

0 . 40000E+01 



-2 

2 

0 . 40000E+01 



-2 

3 

0 . 40000E+01 



-2 

4 

0 . 40000E+01 



-2 

5 

0 . 40000E+01 



-2 

6 

0 . 40000E+01 



-1 

2 

8 0 

0 

6 

-2 

1 

0 . 40000E+01 



-2 

2 

0 . 40000E+01 



-2 

3 

0 . 40000E+01 



-2 

4 

0 . 40000E+01 



-2 





-3 






the normalized displacements u x ld and u y ld at all nodal points 

100CLCASE1 0 . 12052E-01 Fiber Microbuckling 3 1000 


-4 

DISPLACE 3 

1 0 


-5 

UX 1 

2 10 

0 

-5 

UY 1 

2 2 0 

0 

-5 

U 1 

2 0 0 

1ALL 

-1 

1 0 . 14712E+01 

-0 . 32663E+01 


-1 

2 0 . 80530E+00 

-0 . 32449E+01 


-1 

3 0 . 14974E+00 

-0 . 31895E+01 


-1 

4-0 . 31371E+00 

-0 . 31129E+01 


-1 




-3 





the normalized stresses o x Jx y , o yy /x y , o xy /x y and q, x /x v , the normalized effective shear stress 
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T e /t y , the strains y s and y p , the fiber misalignment cp and others [20] 


100CLCASE1 0.12052E 

-01 


Fiber Microbuckling 3 

1000 


-4 

STRSSTRN 

12 

3 

0 





-5 

SIG_XX 

1 

4 

1 

1 0 




-5 

SIG_YY 

1 

4 

2 

2 0 




-5 

SIG_XY 

1 

4 

1 

2 0 




-5 

SIG_YX 

1 

4 

2 

1 0 




-5 

TAU_E 

1 

1 

2 

1 0 




-5 

E_LP 

1 

1 

2 

1 0 




-5 

E_LM 

1 

1 

2 

1 0 




-5 

E_T 

1 

1 

2 

1 0 




-5 

GAMA_S 

1 

1 

2 

1 0 




-5 

KAPPA 

1 

1 

2 

1 0 




-5 

GAMA_P 

1 

1 

2 

1 0 




-5 

PHI 

1 

1 

2 

1 0 




-1 

1 8 

0 

0 

6 





-2 

1 -0.12328E+02 

0.20839E+00 

-0 . 11766E+00 

-0 . 13105E+00 

0.14941E+01 

-0 . 30784E-02 

-2 

1 -0 . 25050E- 

-02 

0 . 17207E-03 

0 . 20919E-01 

-0 . 57344E-03 

0 . 10263E-01 

0 . 7 6157E+01 

-2 

2 -0 . 11658E+02 

0.44093E+00 

-0 . 74204E+00 

-0 . 73252E+00 

0.74208E+00 

-0 . 28156E-02 

-2 

2 -0 . 25441E- 

-02 

0 . 18100E-02 

0 . 60815E-02 

-0 . 27141E-03 

0 . 13715E-02 

0 . 66904E+01 

-2 

3 -0.13010E+02 

0 . 17288E-01 

-0 . 74737E+00 

-0 . 74508E+00 

0.79505E+00 

-0 . 30924E-02 

-2 

3 -0 . 28856E- 

-02 

0 . 63897E-04 

0 . 72372E-02 

-0 . 20685E-03 

0 . 15640E-02 

0 . 67907E+01 

-2 

4 -0 . 11966E+02 

0.32614E+00 

-0 . 32097E+00 

-0 . 30034E+00 

0.12303E+01 

-0 . 29931E-02 

-2 

4 -0 . 24544E- 

-02 

0 . 11352E-02 

0 . 14406E-01 

-0 . 53866E-03 

0 . 57257E-02 

0 . 72148E+01 


-2 ... 
-3 


During the analysis, data for three boundaries are output to the output file post.d The 
file post. d contains contains seven columns. Column 1 lists the load step number, column 2 
lists load increase (+1) or decrease (-1), columns 3 and 4 list the normalized displacements ujd 
and Uy/d for the selected node on boundary 3, columns 5 and 6 list the net loads on boundary 
sections 1 and 2 as shown in Figure 9. Column 7 lists global loading parameter A for each load 
step, where A is a measure of how close the simulation is to fulfilling the requested stress loading 
discussed in section 3.6. 


0 

1 

2 

3 

4 


1 0 . 00000E+00 
1 -0.24101E-01 
1 -0 . 37684E-01 
1 -0 . 52874E-01 
1 -0 . 69664E-01 


0 . 00000E+00 
-0 . 13840E-01 
-0 . 22066E-01 
-0 . 31723E-01 
-0 . 43129E-01 


0 . 00000E+00 
0 . 10106E+03 
0 . 15836E+03 
0 . 22269E+03 
0 . 29416E+03 


0 . 00000E+00 
0 . 10106E+03 
0 . 15836E+03 
0 . 22269E+03 
0 . 29416E+03 


0 . 00000E+00 
0 . 96809E-03 
0 . 15167E-02 
0 . 21325E-02 
0 . 28163E-02 
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996 

1 

-0 . 49693E+00 

0 . 21751E+00 

0.19896E+04 

0 . 19896E+04 

0 . 18922E-01 

997 

1 

-0 . 49626E+00 

0 . 21723E+00 

0 . 19925E+04 

0 . 19925E+04 

0 . 18949E-01 

998 

1 

-0 . 49530E+00 

0 . 21693E+00 

0 . 19966E+04 

0 . 19966E+04 

0 . 18988E-01 

999 

1 

-0 . 49464E+00 

0 . 21683E+00 

0 . 19995E+04 

0 . 19995E+04 

0 . 19015E-01 

1000 

1 

-0 . 49373E+00 

0 . 21680E+00 

0 . 20036E+04 

0 . 20036E+04 

0 . 19053E-01 


6.2. Result Visualization 

The output file flash . out is formatted to be used with the commercial post-processing 
software FEMGV 3 . A translation routine (make_post . f ) was developed, which reads 
flash, out and creates output data of the model (flash_pat . out) and displacements 
(flash . dis) in neutral file format to be used with PATRAN 4 . The routine can also create output 
of the model and result data (flash.dat) to be used with Tecplot 5 . The routine can easily be 
modified to create an output file in a particular format that is compatible with a graphics package 
of the user’s choice. 


6.2.1 Examples of Deformation Plots 

The deformed finite element meshes for the unit cell under axial compression are shown in 
Figure 19. At the beginning of the loading a symmetric deformation with respect to the two axis 
of symmetry can be observed as shown in Figure 19a. At the edge where the axial compression 
stress is applied a non-uniform displacement distribution can be observed which forms an axial 
contraction. The increased displacement is caused by a reduced stiffness which is the result of 
modeling the resin pocket and z-pin as a void. At the end of the analysis after passing the failure 
point and the formation of kink bands the deformation grows to be unsymmetrical and the unit 
cell becomes increasingly distorted as shown in Figure 19b. 

The deformed mesh for axial and transverse compression is shown in Figure 20. The 
transverse compression load was 2% of the compression load applied axially as shown in 
Figure 12. Analysis without constraints as shown in Figure 20 suggested that constraining points 
A and B shown in Figure 12, was not required for small transverse compression loads. Compared 
to the deformation shown in Figure 19a for the pure compression case the axial displacement 
shown in Figure 20 is reduced at the top and bottom edges where the transverse compression is 
applied and the void remains open. 

The deformed mesh for axial compression and transverse tension are shown in Figure 21. 
As for the previous case the transverse load was chosen to be 2% of the compression load applied 
axially as shown in Figure 13. Compared to the deformation for the pure compression case as 
shown in Figure 19a the axial displacement shown in Figure 21 has increased at the top and 
bottom edges where the transverse tension stress is applied; also the void opening has increased. 


3 http://www.femsys.co.uk/ 

4 http://www.mscsoftware.com/products/quick_prod.cfm 

5 http://www.tecplot.com/ 


14 



6.2.2 Examples of Contour Plots 

As an example the contour plots of the fiber misalignment (f> , the effective shear stress x e and the 
strains y s and y p , are shown in Figures 24-27 for the unit cell subjected to axial and transverse 
compression. The contours correspond to the end of the analysis after passing the failure point 
where the severe color gradients indicate the formation and location of kink bands. 


6.2.3 Examples of Stress-Displacement Plots to Obtain Compression Strength 

As an example typical stress versus displacements plots are shown in Figures 28 and 29. Typically 
the post . d file is read into a graphing program such as KaleidaGraph™ or Excel™. Here the 
stress is obtained by multiplying the global loading parameter A from column (#7) of the output 
file post . d with the matrix yield stress r y from the Ramberg-Osgood strain hardening law that is 
used as input and the requested stress loading (-1000) described in section 3.6. The compression 
stress o xx = A r y 1000 is plotted versus the normalized displacement in x-direction (u jd) for the 
selected node on boundary 3 as shown in Figure 9. The normalized displacement (ujd) is listed in 
column (#3) of the output file post.d. As shown in Figures 28 and 29, the plot typically indicates a 
clear first maximum. The exact maximum value, which is the predicted compression strength, is best 
found in the data rather than the plot. 
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LIST OF TABLES 


Table 1. 

Dimensions of Unit Cell 


Graphite/Epoxy UD Prepreg with 2% Z-Pins [13] 


from 

normalized with d ( 

D z 

0.28 mm 


D’z 

0.3 mm 

42.8 

H z 

1.75 mm 

250 

L z 

1.75 mm 

250 

C 

1.29 mm 

184 

d 

7 pm 

1 
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Table 2. 

Material Properties for Initial Preliminary Analysis 


Graphite/Epoxy Prepreg (A) [13] 


from 

normalized with r xy 

E u 

117 GPa 

1083 

E 11 (tension) 

9.0 GPa 

83.33 

E 02 (compression) 

9.5 GPa 

87.96 

G n 

4.8 GPa 

44.44 

Gf 

22 GPa 

203.7 

Xy 

108 MPa 

1 

d 

7 pm 

- 

V f 

0.55 

- 

a 

3/7 

- 

n 

3 

- 

Generic Graphite/Epoxy Prepreg (B) [13] 

normalized with x xy 

E u 

- 

1200 

E 12 (tension) 

- 

74 

E 12 (compression) 

- 

87 

G n 

- 

41 


- 

175 

Xy 

- 

- 

d 

7 pm 

- 

Vf 

0.63 

- 

a 

0.429 

- 

n 

3 

- 
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Table 3. 


Carbon/Epoxy Material Properties 


Graphite/Epoxy UD Prepreg 


from 

normalized with r xy 

E u 

161 GPa 

4248 

E 01 (tension) 

11.4 GPa 

301 

E 22 (compression) 

12.8 GPa 

338 

G n 

5.17 GPa 

136 

Gt 

22 GPa 

580 

Ty 

37 MPa 

1 

d 

5.1 pm 

- 

Vf 

0.59 

- 

a 

0.00923 

- 

n 

8.54 

- 
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Table 4. 

Dimensions of Unit Cells 


Graphite/Epoxy UD Prepreg with 2% large diameter Z-Pins 


from 

normalized with d 

D z 

0.508 mm 

- 

D’z 

0.528 mm 

103.53 

H z 

3.175 mm 

622.55 


3.175 mm 

622.55 

C 

2.1844 mm 

428.31 

Graphite/Epoxy UD Prepreg with 4% small diameter Z-Pins 


from 

normalized with d 

D z 

0.28 mm 

- 

D’ z 

0.3 mm 

58.8 

H z 

1 .2446 mm 

244 

L z 

1 .2446 mm 

244 

C 

0.868 mm 

170.2 

Graphite/Epoxy UD Prepreg with 2% small diameter Z-Pins 


from 

normalized with d 

D z 

0.28 mm 

- 

D’ z 

0.3 mm 

58.8 

H z 

1.7526 mm 

343.65 

L z 

1.7526 mm 

343.65 

C 

0.868 mm 

170.2 
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Table 5. 


Load Cases for Laminate Analysis and Results 


external load N x =- 1000 lbs/in 


On, psi 

O 22 , psi 

T 12 , psi 

[0/90] s 

-1898 

-29.5 (-2% a n ) 

0 

[0/±45] s 

-2304 

+48.9 (~2% On) 

0 

[0/45/-45/90] s 

-2597 

-3.3 (-0.2% a u ) 

0 

external loadN x =-1000 lbs/in, N xy =-1000 lbs/in 


On, psi 

o 2 2 , psi 

T 12 , psi 

[0/90] s 

-1898 

-29.5 

-1000 (-50% On) 


-2304 

+48.9 (-2% o n ) 

-222 (-10% o n ) 

[0/45/-45/90] s 

-2597 

-3.3 

-275 (-10% o n ) 


Table 6. 


FLASH Input for Load Cases Used for Strength Reduction Analysis 



axial 

compression 

compression/ 

2%transverse 

tension 

compression 
10% shear 

compression 
50% shear 

compression 
2% tension 
10% shear 

On/ Ty 

-1000 

-1000 

-1000 

-1000 

-1000 

022/f 

- 

+20 

- 

- 

+20 

r 12/f 

- 

- 

100 

500 

100 

Hl/ty 

- 

- 

100 

500 

100 
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area of maximum 
waviness on flanks of 



Figure 1: Typical shape of local area around z-pin 


22 


size of resin pocket 



C 

D z pin diameter 
D’ z pin + resin 
H z vertical spacing 
L z horizontal spacing 
d fiber diameter 


Figure 2: Detail of perfect z-pin field with typical dimensions and unit cell 



z-pin 


resin pocket 


Figure 3: Micrograph of individual z-pin with surrounding laminate [13] 


23 


H z /f 




y 

i 




. C/f 
L z /f 



Figure 4: Typical finite element mesh of unit cell with normalized dimensions and 
initial misalignment due to z-pin insertion 
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24 


25 
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1 

0.00 

0 . 00 

2 

125.00 

0 . 00 

3 

250.00 

0 . 00 

4 

0.00 

103.60 

5 

125.00 

103.60 

6 

250.00 

103.60 

7 

110.00 

105.00 

8 

140.00 

105.00 

9 

0.00 

124 . 00 

10 

33.00 

124 . 00 

11 

217.00 

124 . 00 

12 

250.00 

124 . 00 

13 

32.00 

125.00 

14 

217.99 

125.00 

26 

250.00 

250 . 00 


Figure 5: Outline of unit cell and void with geometry points 
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25 
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5, 
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0 
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0 . 1 
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8, 

11 
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6 

, 12, 
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.25, 

.25 

r 
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10 

9, 

9, 

15 

, 10 

/ 

13, 16 

, 3, 

10 

, 1, 

1, 

o . 

1, 0 

. 1 

11, 

, 14 

r 

17, 

12 

, 12, 

18, 

3, 

10, 

1, 1 

t 

10, 

10 

17, 

, 20 

r 

22, 

18 

, 18, 

23, 

8, 

10, 

4, 4 

r 

10, 

10 

21, 

, 21 

r 

24, 

22 

, 22, 

25, 

8, 

10, 

10, 

10 

, 0. 

1, 0.1 


Figure 6: Outline of unit cell and void with geometry points and mesh sets 
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(b). Contour map of fiber alignment angle [13] 
Figure 7: Measurement of fiber misalignment angle 
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26 


1 

0 . 

. 00 

0 . 

. 00 


2 

125 . 

. 00 

0 . 

. 00 


3 

250 . 

. 00 

0 . 

. 00 

^MISALIGNMENT TYPE 

4 

0 . 

. 00 

103. 

. 60 

2 

5 

125 . 

. 00 

103. 

. 60 

^MISALIGNMENT 

6 

250 . 

. 00 

103. 

. 60 

10,10,26.0,26.0 

7 

110 . 

. 00 

105. 

. 00 
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250 . 

. 00 

250 . 

. 00 



Figure 8: Outline of unit cell with geometry points and fiber misalignment data points 
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u x =0.0 



*STRESS LOADING 
10, 1 


1, 1, 

4, 1 

r 

i 


-1000, 

0, 

o. 

0, 

0 

4, 4, 

9, 1 

t 

1 


-1000, 

0, 

0, 

0, 

0 

9, 9, 

15, 

1, 

1 


-1000, 

0, 

0, 

0, 

0 

15, 15 

, 21 

r 

1, 

1 

-1000, 

0, 

o. 

0, 

0 

21, 21 

, 24 

r 

1, 

1 

-1000, 

0, 

o. 

0, 

0 

3, 3, 

6, 1 

r 

1 


-1000, 

0, 

o. 

0, 

0 

6, 6, 

12, 

i. 

1 


-1000, 

0, 

o. 

0, 

0 

12, 12 

, 18 

r 

1, 

1 

-1000, 

0, 

o. 

0, 

0 

18, 18 

, 23 

r 

1, 

1 

-1000, 

0, 

o. 

0, 

0 

23, 23 

, 26 

r 

1, 

1 

-1000, 

0, 

o. 

0, 

0 

Figure 9 

: Finite 

elemei 


*TYPE 1 CONSTRAINTS 
3 


1, 

1, 2, 21, 1, 
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Figure 10: M esh with axial compression load and boundary contionsfor UD Ply without Z-Pin 
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Figure 11: Calculation of ply stresses using classical laminate theory 


31 


U X =0.0 Gyy 



U X =0.0 Gyy 


^STRESS LOADING 
14, 1 

1, 1, 4, 1, 1 
- 1000 , 0 , 0 , 0 , 0 

4, 4, 9, 1, 1 
- 1000 , 0 , 0 , 0 , 0 

9, 9, 15, 1, 1 

- 1000 , 0 , 0 , 0 , 0 


*TYPE 1 CONSTRAINTS 
3 


1, 1, 2, 21, 

1, 

1, 

0 


9, 9, 15, 4, 

4, 

2, 

0 


24, 24, 25, 

21, 

1, 

i. 

0 

5, 8, 11, 1 

, 17 

r 

2, 

0 

16, 19, 22, 

17, 

i. 

2, 

0 


23, 23, 26, 1, 1 
- 1000 , 0 , 0 , 0 , 0 
1 , 1 , 2 , 1 , 1 
0 , - 20 , 0 , 0 , 0 

2, 2, 3, 1, 1 

0 , - 20 , 0 , 0 , 0 

24.24.25, 1, 1 

0 , - 20 , 0 , 0 , 0 

25.25.26, 1, 1 

0 , - 20 , 0 , 0 , 0 

Figure 12: Finite element model subjected to axial and transverse compression loading 


32 


Oxx 

Uy=0.0 


y 

1 


Ux=0.0 CTyy 



Ux=0.0 (5yy 

^STRESS LOADING 
14, 1 


1, 1, 

4, 

1, 

1 


*TYPE 

1 CONSTRAINTS 


-1000, 

0, 

0, 

0, 

0 

3 





4, 4, 

9, 

1, 

1 


1, 

1, 2, 21, 

1, 

1, 

0 

-1000, 

0, 

0, 

0, 

0 

9, 

9, 15, 4, 

4, 

2, 

0 

9, 9, 

15, 

1, 

1 


24, 

24, 25, 

21, 

1, 

1, 

-1000, 

0, 

0, 

0, 

0 






15, 15 

, 2 

1, 

1, 

1 






-1000, 

0, 

0, 

0, 

0 







23, 

23, 

VO 

CM 

1, 1 

-1000, 

O 

O 

o 

o 

1, 

1, 2 

, 1, 

1 

0, 

o 

C\j 

o 

o 

0 

2, 

2, 3 

, 1, 

1 

0, 

20, 

o 

o 

0 


24.24.25, 1, 1 
0 , 20 , 0 , 0 , 0 

25.25.26, 1, 1 
0 , 20 , 0 , 0 , 0 


Oxx 


0 


Figure 13: Finite element model subjected to axial compression and transverse tension loading 
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Figure 14: Finite element model subjected to axial compression and shear loading 
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Figure 15: Finite element model subjected to compression, transverse tension and shear loading 
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Figure 16:Finite element mesh of unit cell for small z-pin with 2% areal density 
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Figure 17:Finite element mesh of unit cell for small z-pin with 4% areal density 
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Figure 18:F i nite element mesh of unit cel I for large z-pin with 2% areal density 
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undeformed 

outline 


Figure 19: Analysis of unit cell of UD Ply with Z-Pin subjected to axial compression 
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outline 



Figure 20: Deformed finite element model of unit cell subjected to 
axial and transverse compression loading 



Figure 21: Deformed finite e 


of unit cell subjected to 


axial compression and transverse tension loading 
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undeformed 

outline 



Figure 22: Deformed finite element model of unit cell subjected to 
axial compression and shear loading 



Figure 23: Deformed finite element model of unit cell subjected to 
axial compression, transverse tension and shear loading 
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Figure 24: Contour plot of fiber misalignment angle 0 
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Figure 25: Contour plot of shear stress T e 
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Figure 26: Contour plot of shear strain y p 
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Figure 27: Contour plot of shear strain y s 
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